% Copyright 2020 Makani Technologies LLC
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%      http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

function power_curve(mean_wind_speed)
  params = make_params();

  %% figure('name', 'power curve');
  %% hold on;
  %% wind_speeds = linspace(0, 25, 1000);

  %% min_kite_airspeed = 30;
  %% max_tension = 250e3;
  %% kite_area = 32.9;
  %% max_grid_power = 600e3;
  %% max_zeta = calc_max_zeta(max_tension, kite_area, max_grid_power, params)

  %% m600_capacity_factor = ...
  %%     calc_capacity_factor(mean_wind_speed, min_kite_airspeed, max_tension, ...
  %%                          max_zeta, kite_area, max_grid_power, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, min_kite_airspeed, ...
  %%                                 max_tension, max_zeta, kite_area, ...
  %%                                 max_grid_power, params);

  %% calc_tether_drag_coeff(max_tension, kite_area, max_grid_power, params)

  %% critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
  %%                                                  max_tension, ...
  %%                                                  max_zeta, kite_area, ...
  %%                                                  max_grid_power, ...
  %%                                                  params);
  %% power_curve = calc_simple_power_curve(wind_speeds, min_kite_airspeed, ...
  %%                                       max_tension, max_zeta, kite_area, ...
  %%                                       max_grid_power, params);
  %% plot(wind_speeds, power_curve / 1e3);
  %% %plot(wind_speeds, weibull(wind_speeds, mean_wind_speed, params.weibull_k) .* power_curve / 1e2);

  %% min_kite_airspeed = 30;
  %% max_tension = 0.75 * 250e3;
  %% kite_area = 32.9;
  %% max_grid_power = 600e3;
  %% max_zeta = calc_max_zeta(max_tension, kite_area, max_grid_power, params)

  %% reduced_tension_capacity_factor = ...
  %%     calc_capacity_factor(mean_wind_speed, min_kite_airspeed, max_tension, ...
  %%                          max_zeta, kite_area, max_grid_power, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, min_kite_airspeed, ...
  %%                                 max_tension, max_zeta, kite_area, ...
  %%                                 max_grid_power, params);
  %% calc_tether_drag_coeff(max_tension, kite_area, max_grid_power, params)


  %% critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
  %%                                                  max_tension, ...
  %%                                                  max_zeta, kite_area, ...
  %%                                                  max_grid_power, ...
  %%                                                  params);
  %% power_curve = calc_simple_power_curve(wind_speeds, min_kite_airspeed, ...
  %%                                       max_tension, max_zeta, kite_area, ...
  %%                                       max_grid_power, params);
  %% plot(wind_speeds, power_curve / 1e3);
  %% %plot(wind_speeds, weibull(wind_speeds, mean_wind_speed, params.weibull_k) .* power_curve / 1e2);

  %% min_kite_airspeed = 30;
  %% max_tension = 0.5 * 250e3;
  %% kite_area = 1.2 * 32.9;
  %% max_grid_power = 600e3;
  %% max_zeta = calc_max_zeta(max_tension, kite_area, max_grid_power, params)

  %% half_tension_capacity_factor = ...
  %%     calc_capacity_factor(mean_wind_speed, min_kite_airspeed, max_tension, ...
  %%                          max_zeta, kite_area, max_grid_power, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, min_kite_airspeed, ...
  %%                                 max_tension, max_zeta, kite_area, ...
  %%                                 max_grid_power, params);


  %% calc_tether_drag_coeff(max_tension, kite_area, max_grid_power, params)
  %% critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
  %%                                                  max_tension, ...
  %%                                                  max_zeta, kite_area, ...
  %%                                                  max_grid_power, ...
  %%                                                  params);
  %% power_curve = calc_simple_power_curve(wind_speeds, min_kite_airspeed, ...
  %%                                       max_tension, max_zeta, kite_area, ...
  %%                                       max_grid_power, params);
  %% plot(wind_speeds, power_curve / 1e3);
  %plot(wind_speeds, weibull(wind_speeds, mean_wind_speed, params.weibull_k) .* power_curve / 1e2);

  %% plot_power_curve(wind_speeds, 30, 250e3, 42, 32.9, 600e3, params);
  %% calc_capacity_factor(mean_wind_speed, 30, 250e3, 42, 32.9, 600e3, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, 30, 250e3, 42, 32.9, 600e3, params)

  %% plot_power_curve(wind_speeds, 30, 0.75 * 250e3, 42, 32.9, 600e3, params)
  %% calc_capacity_factor(mean_wind_speed, 30, 0.75 * 250e3, 42, 32.9, 600e3, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, 30, 0.75 * 250e3, 42, 32.9, 600e3, params)

  %% plot_power_curve(wind_speeds, 30, 0.5 * 250e3, 42, 32.9 * 1.2, 600e3, params)
  %% calc_capacity_factor(mean_wind_speed, 30, 0.5 * 250e3, 42, 32.9 * 1.2, 600e3, params)
  %% calc_analytical_capacity_factor(mean_wind_speed, 30, 0.5 * 250e3, 42, 32.9 * 1.2, 600e3, params)

  %% ylabel('Power [kW]');
  %% ylim([0, 700]);
  %% yyaxis right;

  %% plot(wind_speeds, weibull(wind_speeds, mean_wind_speed, params.weibull_k));
  %% ylim([0, 0.15]);
  %% ylabel('Probability density [(m/s)^{-1}]')

  %% grid on;
  %% xlabel('Wind speed [m/s]');
  %% title('Power curve comparison');

  %% legend(sprintf('M600 (C = %0.3f)', m600_capacity_factor), ...
  %%        sprintf('3/4 tension (C = %0.3f)', reduced_tension_capacity_factor), ...
  %%        sprintf('1/2 tension, 1.2x area (C = %0.3f)', half_tension_capacity_factor), ...
  %%        'Wind distribution, Class III', ...
  %%        'Location', 'east');

  %% return;

  % Tunable parameters.
  max_grid_power = 900e3;                     % [W]
  nominal_tension = 150e3;                    % [N]
  nominal_tether_diameter = ...               % [m]
      calc_tether_diameter(nominal_tension, max_grid_power, params);
  min_kite_airspeed = 25;                     % [m/s]
  max_tensions = linspace(50e3, 500e3, 100);  % [N]
  kite_areas = linspace(10, 100, 50);         % [m^2]
  [max_tension_grid, kite_area_grid] = meshgrid(max_tensions, kite_areas);
  max_zeta_grid = calc_max_zeta(max_tension_grid, kite_area_grid, ...
                                max_grid_power, params);
  system_drag_coeff_grid = calc_system_drag_coeff(max_tension_grid, ...
                                                  kite_area_grid, ...
                                                  max_grid_power, params);
  tether_drag_coeff_grid = calc_tether_drag_coeff(max_tension_grid, ...
                                                  kite_area_grid, ...
                                                  max_grid_power, params);

  capacity_factor = zeros(length(kite_areas), length(max_tensions));
  for i = 1:length(kite_areas)
    for j = 1:length(max_tensions)
      capacity_factor(i, j) = calc_capacity_factor(mean_wind_speed, ...
                                                   min_kite_airspeed, ...
                                                   max_tension_grid(i, j), ...
                                                   max_zeta_grid(i, j), ...
                                                   kite_area_grid(i, j), ...
                                                   max_grid_power, ...
                                                   params);
    end
  end
  analytical_capacity_factor = calc_analytical_capacity_factor( ...
      mean_wind_speed, min_kite_airspeed, max_tension_grid, max_zeta_grid, ...
      kite_area_grid, max_grid_power, params);

  figure('name', '\zeta');
  max_zeta_cplot = contour(max_tensions / 1e3, kite_areas, max_zeta_grid);
  clabel(max_zeta_cplot);

  figure('name', 'Tether drag coefficient');
  tether_drag_coeff_cplot = contour(max_tensions / 1e3, kite_areas, ...
                                    tether_drag_coeff_grid);
  clabel(tether_drag_coeff_cplot);

  figure('name', 'System drag coefficient');
  system_drag_coeff_cplot = contour(max_tensions / 1e3, kite_areas, ...
                                    system_drag_coeff_grid);
  clabel(system_drag_coeff_cplot);

  figure('name', '\zeta');
  max_zeta_cplot = contour(max_tensions / 1e3, kite_areas, max_zeta_grid);
  clabel(max_zeta_cplot);

  figure('name', 'Capacity factor');
  ax1 = axes('Position', [0, 0, 1, 1], 'Visible', 'off');
  params_str = sprintf( ...
      ['{\\bfWind parameters}:\n', ...
       'mean wind speed: %0.1f m/s\n', ...
       '\\rho_{air}: %0.3f kg/m^3\n', ...
       'Weibull k: %0.02f\n', ...
       'shear exponent: 0.0\n\n', ...
       '{\\bfPower parameters}:\n', ...
       'max grid power: %0.0f kW\n', ...
       '\\eta_{air-grid}: %0.2f\n\n', ...
       '{\\bfKite parameters}:\n', ...
       'C_L: %0.2f\n', ...
       'C_D_p: %0.3f\n', ...
       'C_D_{ind}: %0.3f\n\n', ...
       '{\\bfTether parameters}:\n', ...
       'length: %0.0f m\n', ...
       'diameter @ %0.0f kN: %0.03f m\n', ...
       'sectional C_D: %0.1f\n\n', ...
       '{\\bfPath parameters}:\n', ...
       'elevation: %0.0f deg\n', ...
       'min airspeed: %0.0f m/s\n', ...
       'gravity power loss: %0.0f kW'], ...
       mean_wind_speed, ...
       params.air_density, ...
       params.weibull_k, ...
       max_grid_power / 1e3, ...
       params.air_to_grid_efficiency, ...
       params.lift_coeff, ...
       params.parasitic_drag_coeff, ...
       params.induced_drag_coeff, ...
       params.tether_length, ...
       nominal_tension / 1e3, ...
       nominal_tether_diameter, ...
       params.tether_sectional_drag_coeff, ...
       params.flight_path_elevation * 180/pi, ...
       min_kite_airspeed, ...
       params.gravity_power_loss / 1e3);
  text(0.72, 0.9, params_str, 'FontSize', 8, 'VerticalAlignment', 'top');
  ax2 = axes('Position', [0.1, 0.1, 0.6, 0.8]);
  capacity_factor_plot1 = contour(ax2, max_tensions / 1e3, kite_areas, ...
                                  capacity_factor, linspace(0, 1, 21));
  %% hold on;
  %% capacity_factor_plot2 = contour(max_tensions / 1e3, kite_areas, ...
  %%                                 analytical_capacity_factor, ...
  %%                                 linspace(0, 1, 21), '--');
  %% legend('Numerical', 'Analytical');
  clabel(capacity_factor_plot1);
  %% clabel(capacity_factor_plot2);
  test_tension = 150e3;  % [N]
  title(sprintf('Capacity factor, Class III (%0.1f m/s)', mean_wind_speed));
  xlabel('Max. tension [kN]');
  ylabel('Kite area [m^2]');
  grid on;
end

function params = make_params()
  % Physical parameters.
  params.air_density = 1.225;  % [kg/m^3]
  params.weibull_k = 2.0;      % [#] 2.0 is equivalent to Rayleigh distribution.

  % Power system parameters.
  params.air_to_grid_efficiency = 0.67;  % [#]

  % Kite parameters.  The lift and span efficiency numbers assume a equivalent
  % monoplane CL=2 with a knock-down for a biplane configuration
  % (4 m separation) estimated from AVL.
  params.lift_coeff = 2.0;              % [#]
  params.aspect_ratio = 20;             % [#]
  params.span_efficiency = 0.66;        % [#]
  params.parasitic_drag_coeff = 0.055;  % [#]
  params.induced_drag_coeff = params.lift_coeff^2 / ...
      (pi * params.aspect_ratio * params.span_efficiency);

  % Tether parameters.
  params.tether_length = 430;                                   % [m]
  params.tether_sectional_drag_coeff = 0.7;                     % [#]
  params.use_nominal_tether_diameter = false;                   % [Boolean]
  params.tether_diameter_power_coeff = 0.014 / sqrt(600e3);     % [m/sqrt(W)]
  params.tether_diameter_tension_coeff = 0.0052 / sqrt(250e3);  % [m/sqrt(N)]
  params.tether_diameter_padding = 0.0074;                      % [m]
  params.nominal_tether_diameter = ...                          % [m]
      params.tether_diameter_power_coeff * sqrt(600e3) + ...
      params.tether_diameter_tension_coeff * sqrt(250e3) + ...
      params.tether_diameter_padding;

  % Flight path parameters.
  params.flight_path_elevation = 30 * pi/180;  % [rad]
  params.gravity_power_loss = 90e3;            % [W] This should be positive.
end

function tether_diameter = calc_tether_diameter(max_tension, max_grid_power, ...
                                                params)
  if params.use_nominal_tether_diameter
    tether_diameter = params.nominal_tether_diameter;
  else
    % Approximate tether diameter scaling based on construction details at:
    % go/makanitether.
    tether_diameter = ...
        params.tether_diameter_power_coeff * sqrt(max_grid_power) + ...
        params.tether_diameter_tension_coeff * sqrt(max_tension) + ...
        params.tether_diameter_padding;
  end
end

function tether_drag_coeff = calc_tether_drag_coeff(max_tension, kite_area, ...
                                                    max_grid_power, params)
  tether_diameter = calc_tether_diameter(max_tension, max_grid_power, params);
  tether_drag_coeff = 0.25 * params.tether_sectional_drag_coeff * ...
      params.tether_length * tether_diameter ./ kite_area;
end

function system_drag_coeff = calc_system_drag_coeff(max_tension, kite_area, ...
                                                    max_grid_power, params)
  tether_drag_coeff = calc_tether_drag_coeff(max_tension, kite_area, ...
                                             max_grid_power, params);
  system_drag_coeff = params.parasitic_drag_coeff + ...
                      params.induced_drag_coeff + ...
                      tether_drag_coeff;
end

function max_zeta = calc_max_zeta(max_tension, kite_area, max_grid_power, ...
                                  params)
  system_drag_coeff = calc_system_drag_coeff(max_tension, kite_area, ...
                                             max_grid_power, params);
  max_zeta = 4/27 * params.lift_coeff^3 ./ system_drag_coeff.^2;
end

function probability_density = weibull(wind_speed, mean_wind_speed, k)
  probability_density = k / mean_wind_speed * gamma(1 + 1 / k) * ...
      (wind_speed / mean_wind_speed .* gamma(1 + 1 / k)).^(k - 1) .* ...
      exp(-(wind_speed / mean_wind_speed * gamma(1 + 1 / k)).^k);
end

% Numerical method of calculating capacity factor.
function capacity_factor = calc_capacity_factor( ...
    mean_wind_speed, min_kite_airspeed, max_tension, max_zeta, kite_area, ...
    max_grid_power, params)
  wind_speeds = linspace(0, 50, 1000);
  dwind = wind_speeds(2) - wind_speeds(1);
  wind_distribution = weibull(wind_speeds, mean_wind_speed, params.weibull_k);
  simple_power_curve = calc_simple_power_curve(wind_speeds, ...
                                               min_kite_airspeed, ...
                                               max_tension, max_zeta, ...
                                               kite_area, max_grid_power, ...
                                               params);
  capacity_factor = sum(wind_distribution .* simple_power_curve * dwind) / ...
                    max_grid_power;
end

% Analytic method of calculating capacity factor.
function capacity_factor = calc_analytical_capacity_factor( ...
    mean_wind_speed, min_kite_airspeed, max_tension, max_zeta, kite_area, ...
    max_grid_power, params)
  critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
                                                   max_tension, max_zeta, ...
                                                   kite_area, ...
                                                   max_grid_power, params);
  tension_wind_ratio = critical_wind_speeds.tension / mean_wind_speed;
  power_wind_ratio = critical_wind_speeds.power / mean_wind_speed;
  cutin_wind_ratio = critical_wind_speeds.cutin / mean_wind_speed;

  mean_power_zeta = 0.5 * params.air_density * max_zeta .* kite_area * ...
      mean_wind_speed^3 * params.air_to_grid_efficiency * ...
      cos(params.flight_path_elevation)^3 .* ...
      (6/pi * erf(sqrt(pi)/2 * tension_wind_ratio) - ...
      tension_wind_ratio .* exp(-pi/4 * tension_wind_ratio.^2) .* ...
      (6/pi + tension_wind_ratio.^2));
  mean_power_tension = params.air_to_grid_efficiency * ...
      cos(params.flight_path_elevation) * mean_wind_speed * max_tension .* ...
      (erf(sqrt(pi)/2 * power_wind_ratio) - ...
       erf(sqrt(pi)/2 * tension_wind_ratio) + ...
       (2/3 * tension_wind_ratio - power_wind_ratio) .* ...
       exp(-pi/4 * power_wind_ratio.^2) + ...
       1/3 * tension_wind_ratio .* exp(-pi/4 * tension_wind_ratio.^2));
  mean_power_power = max_grid_power * exp(-pi/4 * power_wind_ratio.^2);
  mean_power_grav = params.gravity_power_loss * ...
      (exp(-pi/4 * power_wind_ratio.^2) - exp(-pi/4 * cutin_wind_ratio.^2));

  mean_power = mean_power_zeta + mean_power_tension + mean_power_power + ...
      mean_power_grav;
  capacity_factor = mean_power / max_grid_power;
end

function plot_power_curve(wind_speeds, min_kite_airspeed, max_tension, ...
                          max_zeta, kite_area, max_grid_power, params)
  critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
                                                   max_tension, ...
                                                   max_zeta, kite_area, ...
                                                   max_grid_power, ...
                                                   params);
  power_curve = calc_simple_power_curve(wind_speeds, min_kite_airspeed, ...
                                        max_tension, max_zeta, kite_area, ...
                                        max_grid_power, params);
  plot(wind_speeds, power_curve / 1e3, '-');
end

function plot_capacity_factor()
  figure(2);
  contour_plot = contour(max_tension / 1e3, kite_area, capacity_factor, ...
                         linspace(0, 1, 21));
  clabel(contour_plot);
  title(sprintf(['Capacity factor: %0.0f kW @ %0.1f m/s\n', ...
                 '%0.2f m tether diameter'], ...
                max_grid_power / 1e3, mean_wind, tether_diameter));
  xlabel('Max. tension [kN]');
  ylabel('Kite area [m^2]');
  grid on;
end

function critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
                                                          max_tension, ...
                                                          max_zeta, ...
                                                          kite_area, ...
                                                          max_grid_power, ...
                                                          params)
  system_drag_coeff = sqrt(4/27 * params.lift_coeff^3 ./ max_zeta);
  min_flight_power = 0.5 * params.air_density * min_kite_airspeed.^3 .* ...
      system_drag_coeff .* kite_area;

  cutin_wind_speed0 = system_drag_coeff / params.lift_coeff .* ...
      min_kite_airspeed / cos(params.flight_path_elevation);
  tension_wind_speed = sqrt(max_tension / 3 ./ (0.5 * params.air_density .* ...
                                                max_zeta .* kite_area)) / ...
                       cos(params.flight_path_elevation);
  critical_wind_speeds = struct( ...
      'cutin0', cutin_wind_speed0, ...
      'cutin', cutin_wind_speed0 .* (1 + params.gravity_power_loss ./ ...
                                         (params.air_to_grid_efficiency * ...
                                          min_flight_power)), ...
      'airspeed', 1.5 * cutin_wind_speed0, ...
      'tension', sqrt(max_tension / 3 ./ ...
                      (0.5 * params.air_density .* max_zeta .* kite_area)) / ...
                 cos(params.flight_path_elevation), ...
      'power', (max_grid_power + params.gravity_power_loss) ./ ...
               (params.air_to_grid_efficiency * max_tension * ...
               cos(params.flight_path_elevation)) + 2/3 * tension_wind_speed);
end

% This "simple" power curve accounts for air-to-grid efficiency
% losses, off-downwind losses, and the airspeed, tension, and
% power constraints.  It does not account for other losses such as
% gravity losses.
function power = calc_simple_power_curve(wind_speed, min_kite_airspeed, ...
                                         max_tension, max_zeta, kite_area, ...
                                         max_grid_power, params)
  critical_wind_speeds = calc_critical_wind_speeds(min_kite_airspeed, ...
                                                   max_tension, ...
                                                   max_zeta, kite_area, ...
                                                   max_grid_power, ...
                                                   params);

  system_drag_coeff = sqrt(4/27 * params.lift_coeff^3 / max_zeta);
  min_flight_power = 0.5 * params.air_density * min_kite_airspeed^3 * ...
      system_drag_coeff * kite_area;

  % Limit power according to airspeed.
  airspeed_limit_power = params.air_to_grid_efficiency * min_flight_power * ...
      (wind_speed / critical_wind_speeds.cutin0 - 1) - ...
      params.gravity_power_loss;
  airspeed_flag = critical_wind_speeds.cutin0 <= wind_speed & ...
                  wind_speed < critical_wind_speeds.airspeed;

  % Limit power according to maximum power take-out.
  zeta_area_limit_power = params.air_to_grid_efficiency * ...
      0.5 * params.air_density * wind_speed.^3 * max_zeta * kite_area * ...
      cos(params.flight_path_elevation)^3 - params.gravity_power_loss;
  zeta_area_flag = critical_wind_speeds.airspeed <= wind_speed & ...
                   wind_speed < critical_wind_speeds.tension & ...
                   wind_speed < critical_wind_speeds.power;

  % Limit power according to tension.
  tension_limit_power = params.air_to_grid_efficiency * max_tension * ...
      cos(params.flight_path_elevation) * ...
      (wind_speed - 2/3 * critical_wind_speeds.tension) - ...
      params.gravity_power_loss;
  tension_flag = critical_wind_speeds.tension <= wind_speed & ...
                 wind_speed < critical_wind_speeds.power;

  % Limit power to maximum grid power at high wind speeds.
  power_limit_power = max_grid_power * ones(size(wind_speed));
  power_flag = wind_speed >= critical_wind_speeds.power;

  power = max(airspeed_limit_power .* airspeed_flag + ...
              zeta_area_limit_power .* zeta_area_flag + ...
              tension_limit_power .* tension_flag + ...
              max_grid_power .* power_flag, 0);
end
